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Modeling of random bimodal structures of 
composites (application to solid propellants): 
I. Simulation of random packs 

V.A. Buryachenko 1 2 Q, T.L. Jackson 23 ), G. Amadio 3 ) 

D ' Abstract: We consider a composite medium, which consists of a homogeneous matrix containing a 

00 ■ statistically homogeneous set of multimodal spherical inclusions. This model is used to represent the 

morphology of heterogeneous solid propellants (HSP) that are widely used in the rocket industry. The 
Lubachevsky-Stillinger algorithm is used to generate morphological models of HSP with large polydisperse 
packs of spherical inclusions. We modify the algorithm by proposing a random shaking procedure that 
leads to the stabilization of a statistical distribution of the simulated structure that is homogeneous, highly 
mixed, and protocol independent (in sense that the statistical parameters estimated do not depend on 
the basic simulation algorithm). Increasing the number of shaking has a twofold effect. First, the system 
becomes more homogeneous and well-mixed. Second, the stochastic fluctuations of statistical parameters 
(such as e.g. radial distribution function, RDF), estimated by averaging of these structures, tend to 
diminish. 

Keywords: Microstructurcs, random packing, spherical particles, inhomogeneous material. 
1. Introduction 

There is a growing recognition within the nano- and micromechanics community that the prediction 
of mechanical properties of heterogeneous materials, such as solid propellants, must take into account 
the various scales (such as the particle sizes) and cannot be treated as a homogeneous material. This 
is due in large part by the use of image analysis and computer-simulation methods on one hand and 
advanced experimental techniques (such as X-ray tomography and electron microscopy) and improved 
materials processing (prescribed structure controlled by processing) on the other. The prediction of the 
behavior of composite materials in terms of the mechanical properties of constituents and their microstruc- 
ture is a central problem of micromechanics, while the quantitative description of the microtopology of 
heterogeneous media is crucial in the prediction of overall mechanical and physical properties of these 
materials. After many years of comprehensive study by direct measurements and empirical relations that 
is extremely laborious, the structure of microinhomogeneous materials is not completely understood. In- 
fluence of microstructure on mechanical properties of microinhomogeneous media drastically increases for 
highly packed particulate media such as, e.g., concentrated suspensions (see, e.g., Chong, Christiansen 
and Baer, 1971), concrete (see for references Garboczi and Bentz, 1997; Amparano, Xib and Roh, 2000), 
dental restorative composite materials (see Lingois and Berglund, 2002; Wissler et ai, 2003; Tanimoto et 
ai, 2004), and solid propellants. 

Heterogeneous solid propellants (HSPs) commonly used in aerospace propulsion are likely to play 
an important role for as long as rockets are built. Fundamentally, they consist of multimodal spherical 
particles of solid oxidizer with extremely high volume fractions (>0.9), typically ammonium pcrchlorate, 
dispersed in a rubbery fuel binder fuel such as hydroxi-terminated-polybutadicnc (HTPB) or polybu- 
tadiene acrylonitrile (PBAN). Designers of rockets are concerned with a number of propellant-related 
issues, including the burning rate, the thermal and mechanical properties of the propellants, and, for 
metallized propellants, the behavior of the metal particles at the surface, including agglomeration, (see 
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Massa, Jackson and Short, 2003; Webb and Davis, 2006; Matous and Geubelle, 2006; Matous et al, 2007; 
Maggi et al, 2008; Xu, Aravas and Sofronis, 2008). 

The prediction of effective mechanical properties of polymer-bonded explosives provides a new chal- 
lenge. These particulate composites contain extremely high volume fractions of explosive particles (>0.9) 
in a polymer binder typically exhibiting viscoelastic behavior (see, e.g., Banerjee and Adams, 2004). In 
explosives, particle sizes range from a few /urn for oxidizer particles in the "dirty binder" to 10s of /jm for 
the metal flakes and 100s of /xm for the largest oxidizer crystals. The oxidizer (HMX) particle shapes are 
more polyhedral than spherical (see Stafford and Jackson, 2010). Large dynamic deformation during pro- 
cessing and shock loading can fracture the microstructure and microcracks can create additional internal 
surfaces or porosity that greatly influences explosive initiation sensitivity (see e.g. Baer et al., 2007). 

In many important cases of practical interest (e.g. solid propellants) , the shape of the heterogeneities 
can be assumed as a spherical. Computer simulation of topologically disordered structures by the random 
packing of hard spherical particles in 2-D and 3-D cases has a long history originated in the theory of 
liquids. This problem is closely connected with the known fundamental problem of statistical physics 
- description of the behavior of the particle system with interaction potential of hard spheres (see e.g 
Binder and Heermann, 1997). Random packing of spheres has been studied very extensively due to their 
technological importance, and their capability of modeling the simple liquid, concentrated suspensions, 
amorphous and powder materials. 

Computer simulation of packing problems can be classified into three groups of methods: molecular 
dynamic, Monte Carlo, and dense random packing. Much progress in the theory of the dense random 
packing was achieved by the use of two kinds of methods: the sequential generation models and the 
collective rearrangement models. In the sequential model with the fixed radii of particles proposed by 
Bennet (1972) [see also the modified algorithms Boudreaux and Gregor, 1977; Lu, Ti, and Ishizaki, 1994; 
Kansal, Truskett and Torquato (2000)], so called cluster growth model, a particle being added to the 
surface of particle cluster which grows outwards is placed sequentially to the point closest to the original 
such that the new particle established contact with three existing spheres in the cluster. The key geometric 
procedure by Jerier et al. (2010) consists in placing one sphere in contact with four non-coplanar spheres. 
In the second type of sequential generation model, called the model of "rigid sphere free fall into a virtual 
box" , one particle is dropped vertically each time from the random point onto the surface of an existing 
particle cluster growing upwards (see e.g. Nolan and Kavanagh, 1992; Cesarano, McEuen and Swilcr 
1995; Kondrachuk, Shapovalov and Kartuzov, 1997, Furukawa, Imai and Kurashige, 2000; Kanuparthi, 
et al, 2009). For contrast, dynamic methods assume the reorganization of whole packing due to either 
short or long range interactions between particles. In the collective rearrangement model (CRM, or 
concurrent), N points randomly distributed in a virtual box are assigned both radii and random motion 
vectors. Each sphere is moved until there are no overlaps. Then the radii are increased and the process 
is repeated until any further increase in radii or any displacement of the spheres create overlaps that can 
not be eliminated [the different versions of this method can be found e.g. in Clarke and Willey, 1987; 
Lubachevsky and Stillingcr, 1990; Lubachevsky, Stillinger and Pinson, 1991; Zinchcnko, 1994; Knott, 
Jackson and Buckmaster, 2001; Kochevets et al., 2001; He and Ekere, 2001; see also Oger et at, 1998 
where a detailed description of the advantages and drawbacks of the different algorithms were presented] . 
Classical molecular dynamics can also be considered as a CRM, where particle size is fixed. The CRMs 
satisfactorily reproduce the real packing properties, but require high time-consuming computation. More 
recently, numerical simulations were performed to realize homogeneous and isotropic packing of spheres 
by various methods, for instance, by assuming hypothetical spheres having dual structure whose inner 
diameter defines the true density and the outer one a nominal density (see Jodrey and Tory, 1985). We 
only consider the spherical particles because only a single parameter, the radius, defines a single possible 
type of contact among particles, which can be estimated analytically and easily. 



V.A. Buryachenko, T.L. Jackson, G. Amadio 



3 



The a priori prediction of the maximum packing fraction, c m , for a system of particles is still an 
open question, despite much study. Packing densities for close lattice packing of monomodal particles 
are 7r/\/12 w 0.9069 (triangular) in the case of disks packing into the plane, and 7r/vT8 rj 0.7405 
(fee or hep) in the case of spheres packed into R 3 . However, a well- mixed random packing does not 
self-assemble into one of the periodic arrangements but instead forms a so-called random close-packed 
(RCP) arrangement which is not well defined, protocol dependent (in both the numerical simulation and 
experimental sense, see McGeary, 1961; Berryman 1983; Cheng, Guo and Lay, 2000; Torquato, Truskett 
and Debenetti, 2000; Aste, Saadatfar and Senden, 2005), and can vary in the range 0.59 = 0.64 for 3D 
unimodal spheres. Here the low value corresponds to the loose packs (LRP), while the upper value is a 
characteristic of so-called dense packs (DRP). A shaking process makes it possible to unlock the spheres 
from the jamming configuration and allows them to find the most homogeneous and mixed arrangement. 
So, Nowak et al. (1997) demonstrated by an electromagnetic vibration exciter that under vibrations the 
bead packing evolves from an initial, low-density configuration towards higher density. The measured 
densities are extrapolated to eliminate finite-size effect that enables one to estimate the DRP's density 
c m = 0.6366 ± 0.004 [Stoyan (1998) proposed elegant hypothesis for c m = 2/tt rj 0.63662]. Various 
algorithms have been devised to simulate reordering due to shaking or vibration of dense packing (see, 
e.g., Barker, 1993; Lee et al, 2009) which reduces the volume concentration of the high density jam 
configuration. 

For multimodal systems, small spheres may be placed into the spaces between packed large spheres. 
For the infinitely different sizes of packed spheres, the maximum c m = 0.637 ■ (2 — 0.637) = 0.868. 
The volume fraction ratio £ between large and small spheres is also important, with maximum packing 
obtained at about £ = 0.60 = 0.75 large particles (Shapiro and Probstein 1992; Maggi, et al, 2008, 2010). 
So, for the packs with size ratio A = 6.5 : 1 Maggi et al. (2008) simulated the packing fractions of 0.687 
and 0.731 for the coarse-to-fine volume ratio £ = 0.90 and £ = 0.75, respectively. Trimodal, multimodal, 
and polydisperse systems can reach even higher packing fractions (Zou et al., 2003). However, packing 
configurations containing a wide range of inclusion concentrations have been investigated less and invite 
further investigation. 

The paper is organized as follow. In Section 2, the quantitative descriptors of the bimodal spherical 
dispersion are discussed in order to describe the pattern of particle location as it really exists rather than 
as described by some assumed model. Since random packing structures are strongly dependent on the 
procedure of their generation, in Section 3 we consider a popular algorithm of the collective rearrangement 
model (CRM) and its combination with the shaking procedure (shaking collective rearrangement model, 
SCRM) adapted for obtaining the most homogeneous configurations which are protocol independent. In 
Section 4 we recall the basic concepts defining some classical method of micromechanics and present ex- 
plicit formula for both effective clastic moduli and stress concentrator factor for composites with bimodal 
distribution of spherical particles. In Section 5 one estimates the RDFs of configurations generated by the 
CRM and SCRM for unimodal and bimodal distributions of spherical particles. Diminishing of stochastic 
fluctuations of the RDFs are reached by using of parallel computing. Comparison of these simulations for 
unimodal packs with the real scanned data obtained by X-ray tomography demonstrates an advantage of 
the SCRM with respect to CRM. One performs a detailed analysis of local maxima and minima of RDFs 
for bimodal simulated packs. 

Estimated RDFs are used in an accompanying paper by Buryachenko (2012) for evaluation of both 
the effective elastic moduli and stress concentrator factors of local stresses in linear elastic composites, 
which consist of a homogeneous matrix containing a statistically homogeneous set of bimodal spherical 
inclusions. 

2 Statistical description of random structure particulate composites 
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Let a linear elastic infinite body R d contains an open bounded domain w C R d (window of obser- 
vation) with a boundary T and with an indicator function W and space dimensionality d (d = 2 and 
d = 3 for 2-D and 3-D problems, respectively). The domain w contains a homogeneous matrix 
and a random finite set X = (wj) (i = 1, . . . ,N(w)) of inclusions w 2 ; with centers Xj and with charac- 
teristic functions Vi(y) equals one if y g u, and zero otherwise and bounded by the spherical surfaces 
T^ 1 ' = {y : |y — Xj| = a' 1 '} and r' 2 ' = {y : |y — Xj| = a^} of two radii a*- 1 ) and a' 2 ', respectively; 
w (fc) = Uv ( k ) (j = 1, 2, ... , AT^, A^ 1 ) + A^ 2 ) = AT, fc = 1, 2). N = N(w) denotes the random number of 
points Xj falling in w; N(w) is the observed number of points in w. 

At first we summarized some basic ideas, notations, and quantities for random point processes such 
as the centers of particles (for more details see Ripley, 1977, 1981: Diggle, 2003; Konig, 1991; Stoyan and 
Stoyan, 1994; Stoyan, 2000, Torquato, 2002). No confusion will arise below in using of terminologies from 
the different scientific subfields. So, in micromechanics one uses the word "field" where statisticians would 
prefer to speak of process although typically there is no time-dependence considered (see, e.g., Stoyan, 
2000). We will consider statistically homogeneous (or stationary) and isotropic ergodic random process 
X, keeping in mind that the stationary, isotropy, and ergodicity can never be tested statistically in their 
full generality. Stationarity means invariance under arbitrary translation, and isotropy means invariant 
under arbitrary rotation. The ergodicity ensured that one sample (one point pattern) is sufficient for 
obtaining statistically secure results, assuming the convergence of results obtained for infinitely expanding 
observation window w. The intensity, n = EX([0, l] d ), of a stationary point process is the mean number 
of points in the unit cube. A standard estimator of the intensity is h = N(w)/w, w = mes w. 

The packed random structure can be characterized by several parameters, such as packing density, 
coordination number, radial distribution function, particle cage, inter-particle spacing, and others. The 
various methods of estimating the effective properties of a composite material use a knowledge of statistical 
geometrical information about the microstructure. In order to incorporate the spatial arrangement of 
components in a micromechanical simulation, it is essential to quantitatively characterize the random 
structure of the composite. The most important factors characterizing the microstructure of composite 
material are the shape, volume fraction, and arrangement (random or regular) of the components that 
permit the calculation of bounds of effective moduli. 

The popular statistical description of microinhomogeneous media is based on expectations of products 
of the characteristic function VW(y), assuming that the role of the matrix is assigned to phase '0'. So, 
improved bounds on a variety of different effective properties have been derived in terms of n-point 
probability density 

sV(y n ) = (V i (y 1 ),...V i (yJ), (1) 

where the angle brackets ((•)) denote an ensemble average; i.e., the probability of simultaneously finding 
of n points in a specified geometrical arrangement y" = y 1; . . . , y„ in one of the phases. In particular, 
two-point probability density is widely used for statistical description of the HSPs (see, e.g., Matous and 
Geubelle, 2006; Maggi et al, 2008; Lee, Gillman and Matous, 2011) 

An alternative related approach of quantitative description of the composite microstructure is based 
on the consideration of the inclusion centers statistically described by the multiparticle probability den- 
sities / Tn (xi, . . . , x m ) that give the probability / m (xi, . . . , x m )cixi . . . dx. m to find a sphere center in the 
vicinities dxi . . . dx. m of the point x m = (xi, . . . , x TO ). The f m are the most basic descriptors that char- 
acterize the structure of many-particle system and have been well-studied in the statistical mechanics of 
liquids (Hansen and McDonald, 1986). In particular, /i = n, where n is the number density of inclusions 
connected with the volume fraction c = nvi, Vi = mesUj. For the statistically homogeneous isotropic me- 
dia in the framework of the simplest "two-point" level, Markov and Willis (1998) demonstrated a simple 
interconnection expressing S^ 1 (y 2 ) as a simple one-tuple integral containing the two-point probability 



V.A. Buryachenko, T.L. Jackson, G. Amadio 



5 



density /2M (r = |xj — x 2 |). Torquato and Stell (1985) have related Sn to multidimensional integrals 
over the infinite set of m-particle densities fi,. ■ ■ f m (rn — > 00). 

The widely used informative function that describes the point distribution is the second- order in- 
tensity function K(r) (called also Ripley's (1977) K function) defined as the number of further points 
expected to be located within a distance r of an arbitrary point divided by the number of points per unit 
area n (Ripley, 1977; Diggle, 1983; see also Pyrz, 2004; Silberschmidt, 2008). Since points lying outside 
the observation window w are nonobscrved, the latter depends strongly on the shape and the size of w. It 
is our aim to simulate a typical realization that also includes interaction to the structure element outside 
of the window while avoiding systematic errors or biases in the estimation procedure. The effect of the 
edge of the domain w becomes increasingly dominant at the dimensional increases. A number of special 
edge corrections are known. A naive way proposed in the minus-sampling method is to consider w* within 
domain w and allow measurements from an object in w* to an object in w. Although the effective sample 
size is then the number of points in w, the method, of course, leads to a big loss of information. A much 
better idea of edge-correction of the estimator for K(r) was suggested by Ripley (1977) 

i j^i 

where N is the number of points x& (fc = 1, . . . N) in the field of observation w with the area w, hj(rij) 
is the indicator function equals 1 if r^ < r and zero otherwise where is a distance between the typical 
fixed points Xj and Xj. The term "typical point" was used by Stoyan and Stoyan (1994) as a synonym 
of "arbitrarily chosen point" , where "arbitrarily" means that a selection scheme is used in which every 
point has the same chance to be selected. The typical point concept is a simplified notion of the Palm 
distribution referring to individual points in a point process (see for details Illian et ai, 2008). Wij is the 
ratio of the circumference contained within w to the whole circumference with radius nj . If the entire 
circle with radius r is placed within the window, Wij = 1 and it is smaller than unity otherwise. For circles 
intersecting the boundary dw, the function Wij compensating for the boundedncss of w is less than one 
and has an explicit formula when the field of observation w is rectangular (see Diggle, 1983; the formula 
and algorithm for 3-D case were proposed by Konig et ai, 1991). The function K(r) (2) is obtained by 
averaging over all inclusions at each value of r. Equation (2) is an approximately unbiased estimator, 
which is free of systematic errors, for sufficiently small r because N/w is a slightly biased estimator for 
n. Illian et al. (2008) considered different improvements for the quality of estimations of K(r) and n. In 
2-D, Diggle (1983) recommended an upper limit of r equal to half the length of the diagonal of a square 
sampling region. 

In a packing simulation only a relatively small number of particles can be considered. To avoid the 
(usually) undesired artificial effects of particles which are not surrounded by neighboring particles (near 
the window of observation boundary) in all directions, one introduces periodic boundary conditions (see 
for details Gajdosfk, and Zeman, Sejnoha, 2006; Steinhauser and Hiermaier, 2009). In such a case, an 
alternative toroidal edge correction (see Ripley, 1979) is often used in which each 2D rectangular box 
w can be regarded as a torus, so that points on opposite edges are considered to be closed (a 3D case 
is analyzed in a similar manner). Then w can be considered to be part of a grid of identical boxes, 
forming a border around the pattern inside w. If a moving particle leaves the central simulation box, 
then one of its image particles enters the central box from the opposite direction. Each of the image 
particles in the neighboring boxes (the numbers of these neighboring boxes equal 8 and 26 in 2D and 
3D case, respectively) moves in exactly the same way. Precisely this method will be explored hereafter 
in this paper for the elimination of the boundary effect. So, we start with a cube w = {[— l,l] d } with 
the center e R d and its periodically located neighbors labeled by the multiplet of integer numbers 
a = (qi, . . . , otd) € Z + where a, = 0, ±1 (i = 1, . . . , d). N random points Xj (i = 1, . . . , N) in the central 
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cube w arc periodically reflected into the neighboring cubes w a . Distances are then measured taking into 
account the minimum image convention according to which the real distance between any two particles 
is given by the shortest distance of any of their images in the surrounding periodic cubes w a (see Ripley, 
1981): |xj — Xj| = min a |x, — x"| = min a |xj — Xj — 2(ai + . . . + a<j)| where x" e w a is a periodic 
mapping of Xj £ w into w a . Toroidal edge correction allows us to analyze all N particles in the central 
box w without using the minus-sampling method because even for the points x, placed near the boundary 
of w all surrounding points are visible. An alternative approach (eliminating the boundary effect of the 
virtual box w) is based on the use of spherical boundary conditions instead of periodic ones (see e.g. Hall, 
1988; Tobochnik and Chapin, 1988). There one simulates hard disks (more exactly a circular cap which 
can be visualized as a contact lens on the surface of an eyeball) on the surface of the ordinary three- 
dimensional sphere and hard spheres on the "surface" of a four-dimensional hypersphere. The advantage 
of this procedure is that there is no preferential direction, and it is impossible to pack particles into 
perfect regular periodic configurations. The considered both toroidal and spherical edge corrections hold 
for simulated structures rather than for the real observation window w obtained by visual cutting out 
from a large sample. In the last case, either the edge correction (2) or minus-sampling method should be 
used. 

For unimodal distributions, the widely used two-point density /2(f) = /2(xi,X2) (r = |xi — X2I) is 
expressed in terms of radial distribution function g(r) (RDF) (called also the pair correlation function, 
see Stoyan and Stoyan, 1994) as 

f 2 (r)=n 2 g(r), (3) 

which is a function of the inter-point distance r. Hence, ng(r)sddr is a probability to find a point whose 
center lies in an infinitesimal spherical shell Sd C R d of inner radius r with a surface Sd (s<j = 27rr, 4wr 2 
for d = 2, 3, respectively) thickness dr and centered in a specific particle. RDF g(r) (3) which plays a role 
similar to the variance in a classical analysis of random variables is defined as the radial distribution of 
the average number of sphere centers per unit area in a spherical shell. The RDF can be estimated from 
second-order intensity function K(r) as (see e.g., Stoyan and Stoyan, 1994; Pyrz, 2004) 

«< r > = dU^^T' (4) 

where u>d is the volume of the unit sphere in R d (o><j = n, 4/3ir for d = 2, 3, respectively). The K- 
function is seldom used in the study of packing, and the RDF g(r) is used instead despite the fact that 
the estimation of g(r) is more delicate due to numerical evaluation of the derivative (4). The RDF is 
related to the derivative of K(r) (4), and is therefore more sensitive to changes in the spatial order 
than is the function K(r). For a completely random point process (i.e. a homogeneous Poisson process) 
K(r) = ujdr d , g(r) = 1. Values of the RDF g(r) larger than 1 indicate that the interposing distances 
around r appear more frequently compared to those in a completely random point process (typically 
there is clustering), whereas values of g(r) smaller than 1 indicate that the corresponding distances are 
rare, that there is mutual inhibition. RDF takes all values from zero to large values at g(r — > a + ) = g(a + ) 
(a + = lim e ->-o ct + e, e > 0), and the second neighbor shell at the large c becomes split exhibiting a 
shoulder. For r — > 00 it tends to 1, and the rate of convergence reflects a degree of randomness of the 
structure considered in a comparison with a completely random Poisson field. The radius ro such that 
g(r) w 1 at r > r determines the range of geometrical disorder where a reference heterogeneity does not 
distinguish individual neighbors and "fills" the surrounding as a continuum. 

Up to now we only considered unimodal distribution of identical heterogeneities. The introduction of 
mark correlation function (see Stoyan and Stoyan, 1994; Gavrikov and Stoyan, 1995; Pyrz, 2004; Rudge, 
Holness and Smith, 2008) may facilitate the characterization of more complex aspect of random structure. 
One attaches to each point with a position x,; of the center set a typical feature which may correspond 
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to the size, shape and type of inclusions served as "point identifiers" . In our case, the "points" are the 
particle locations, while "marks" describe the "particle type" with reference to the particle radii a^- 1 ' and 
a^ 2 \ Then Eq. (2) for the toroidal edge correction (when wtj = 1) can be generalized to the partial K km 
functions 

K km (r) = ^ V V Infra), (5) 

where the indexes i and j correspondent to the particles v^ k ' and t>( m ) fk, m = 1,2), respectively. 
The interpretation of the partial RDF g km is similar to that of g(r). It starts with the probability 
/2|fcm( x i, X2)dxi(ix2 of having one of particle type k and the other of type m. Then (A;, m = 1, 2) 

/ 2 | fcm (x 1)X2 ) = n^n^g km (r), g km (r) = ^-L^^M , (6) 

where n( fc ) is the intensity of particles v^ k \ the mean number fc— points per unit volume, n = n*- 1 ) +n^ 2 \ 
In accordance with the definition, f2\km IS symmetric with respect to k f-> m that implies the symmetry 
gkmfr) = gmkfr) for fc, m = 1, 2 and Vr > + a 1 -™ 1 ' (see e.g., Illian et ai, 2008). Introduction of partial 
descriptors (5) and (6) enables one to describe more accurately the random structures than their mass- 
weighted version of the spatial repartitions proposed by Gallier (2009) for analysis of solid propcllants. 

It is also known other types of statistical descriptors (used, e.g., for the bound estimations of affective 
properties) such as point/g-particle functions, surface-surface correlation functions, nearest-neighbor dis- 
tribution function, linear-path function, two-point cluster function, chord-length distribution function as 
well as the generalized n-point distribution function for the system of identical spheres H n (y m ;y p ~ m ; r q ), 
which is defined to be the correlation associated with finding m points with positions y m on certain sur- 
face within the medium, p — m with positions yP~ m in certain space exterior to the spheres, and q sphere 
centers with positions r q , n = p + q (see for details Pyrz, 1994; Torquato, 2002 and references therein). 
Although higher-order correlation functions (N > 3) are obtained on theoretical grounds, this is not 
a very practicable approach. However, we will estimate the effective elastic properties by the MEFM 
explicitly depending only on the RDFs. Because of this one will consider in details only the RDFs (6). 



3 Packing algorithm 

3.1 Collective rearrangement model 

Collective rearrangement model (CRM) was previously proposed by Lubachevsky, Stillinger, and 
Pinson (1991) for the packs of unimodal spheres and generalized by Knott, Jackson and Buckmaster (2001) 
and Maggi et al. (2008) to polydispcrsc packs (application to heterogeneous solid propcllants). Just for 
completeness we will briefly reproduce the last version of the CRM accompanied by the shaking procedure 
which is adapted to our goal to generate the most homogeneous and highly mixed configurations. 

Let N random points x,; fi = 1, . . . , N) be located in a central box, periodically reflected into neigh- 
boring boxes. Assign initial velocities U; = (uu, . . . , ua) to each random point, with velocity components 
that are independently distributed at random between —1 and +1 and to the uniformly growing inclu- 
sions Vi C {k = 1,2) with the radii a^(t) — a^t. The centers of the inclusions move according to 
the equations 

¥ = Ui (7) 

with a discontinuous change of the vectors at the moment the particle exits through a face of a central 
box as well as during collisions with other inclusions. 
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The collision time is obtained from the condition that the separation distance is the current diameter, 
that for the inclusions vi C v^> and Vj C v^ m ' (k, m = 1, 2) is 

|xj + u;At - Xj - u,-At| = {af ] t + a^k + a { Q k) At + a^At). (8) 

In the case of the collision, the smaller positive root At of Eq. (8) defines the collision time ry. If the 
collision takes place between the inclusions v% and Vj from central and neighboring squares (e.g. identified 
by (-1,1,0)), respectively, then Xj in Eq. (8) should be replaced by Xj + (—1,1,0) with a subsequent 
estimation of a collision time ry . The exit time t[ is estimated as the smallest positive time for exiting 
of the inclusion Vi being considered through one of the sides of the central box. The determination of 
At, = min(Tij , t[ ) allows us to estimate the new inclusion radii do (t+ At* ) and identify either the colliding 
(vi- and Vj*) or exiting («<•) inclusions for which the new velocities and locations are re-estimated. For 
all other inclusions the position vectors are updated according to Eq. (7): x, — s- Xj + VjAt. For these 
inclusions At = At, and velocities Vj remain unchanged. In addition, the collision and exit times are 
corrected ry — > r,j — At, and rf — > rf — At*, respectively. An optimized cell method with analysis of 
only collisions between particles in neighboring cells is used for reduction of the computation time. The 
details of the re-estimation of velocities for the inclusions Wj* and Vj* based on the conservation laws of 
the momentum and energy can be found in Knott Jackson and Buckmaster (2000). 

3.2 Shaking procedure 

A shaking process (called also Monte-Carlo simulation) gives each particle a small random displace- 
ment independent of its neighbors' positions. This makes it possible to unlock the particles from the local 
jamming configuration and allows them to find the most homogeneous and mixed arrangement. Various 
algorithms have been devised to simulate reordering due to shaking or vibration (also called Monte Carlo 
simulation and local-update algorithms, Luijten, 2006) of dense packing (see Barker, 1993) which reduces 
the volume concentration of the high density jam configuration. Packing configurations containing a wide 
range of inclusion concentrations have been investigated less and required some additional consideration. 
In order to describe the algorithm used in this paper, at first one introduces the following definitions. 

For speeding up the calculations, it is necessary to carry out the shaking process only within a local 
shake up window. The random local shaking of a particle Vi C (k = 1,2) is established in a shake 
up window wj^ = {x : |x — x;| < R^a^} whereby the selected particle center x, is randomly moved 
to a new position x^ uniformly distributed in the window with the normalized radius =const. 
(k = 1, 2) depending on the size of spherical particle in the multimodal system. The displacement x^ — x, 
can be chosen as a randomly oriented vector within a sphere with radius Rj^ which permits control 
over the efficiency of the simulation. If the particle does not overlap with any other inclusions the shaking 
is accepted, otherwise the trial shaking is repeated until the number of attempted trial shakings exceed 
some limit . A valid shaking scheme must be ergodic. This means that there is a path in phase space 
from every state to every other state, via a succession of trial moves. 

Another way to speed up the calculations is to check the collision partners Vj (j = j \ , . . . , f" nt ) of 
Vi only in some restricted neighborhood of Xj called testing window w^ m ^ = {x : |x — X j| < r^ m ^} 
(y( fem ) _ _|_ 4. a ( m ) = const., V{ C v^ k \ Vj C t/ m )) which was introduced for unimodal 

particles in 2D case by Buryachenko, et al. (2003). Thus for bimodal packs of particles and we 
have four different sorts of testing windows w[ km ^ (k,m = 1,2). A similar concept of a near-neighbor 
list by Donev, Torquato, and Stillingcr (2005) can be considered as a generalization of testing window 
to the case of nonsphcrical particles. In this model, the neighbors of a particle (near-neighbor list, Af(i)) 
arc considered to be the set of particles whose centers lie within some maximum distance called the 
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testing window w\ of that particle rather than the "geometric neighbors" which can be determined by 
the more computationally expensive Delaunay tessellation (see e.g. Okabe, Boots and Sugihara 1992). 
Donev, Torquato and Stillinger (2005) demonstrated advantage of the near-neighbor list method with 
respect to a traditional cell method (see e.g. Allen and Tildesley, 1987) used for speeding up for the 
neighbor search. One wants to have the testing window w^™ 1 ^ as large as possible so that there is more 
room for the particle u,; C to move without leaving of w[ km \ However, for the larger w[ km \ the 
more neighbors there will be to examine. The optimal balance, as determined by the choice of w[ km * > 
(k, m = 1,2), is chosen numerically at the testing of concrete packs. Introduction of the testing window 
w V* m ) a llows one to improve an efficiency of the shaking simulation if the displacement — is defined 
by possible collision of Vi C with Vj C w[ km \ Namely, let xj — x^ = u'At' (|v'| = 1) be a displacement 
of x.; in I's trying of x^ random shaking where At' = maxj Ai^-, Vj C w^™ 1 ^ and At' is obtained from a 
particular case of Eq. (8) (j € Af(i); 1 = 1,..., N ( s k) ) 

\7c i + vZAt* j -x j \ = aW+a ( - m \ (9) 

where the smaller positive root Atj of Eq. (9) defines the collision time (or, by other words, the amplitude 
of I's shaking) of inclusions Vi and Vj . Found At' with the corresponding moving directions u' allows one to 
find both the random amplitude of shaking uniformly distributed on the interval [0, r] (r = min(At, R^ ), 
At = max; At') with the corresponding direction u' of u,-s moving. Comparison of At > R^J (or <) makes 
it possible to increase (or decrease, respectively) the size i?^ of the shaking window. 

Only the near neighbor set of inclusions Vj (j = j\, . . . , j % n A Tit depends on k,m) which are located 

in the testing window w[ k centered at the point x, (u, C v^ k \ Vj C v^ m \ k,m = 1, 2) are checked for 
overlap, which also reduces the computer time. One determines the optimal size for the shake up window 
R^a^ which provides the minimum average number of trial shaking attempts. This size of shake up 
window provides the fastest stabilization of statistical parameter estimations. However, the values of these 
parameters do not depend on the size of the shaking window and are, furthermore, protocol independent 
(in contrast with the known methods of random close packing simulations, see e.g. Torquato, Truskett 
and Debenetti, 2000). The shaking process passes through all inclusion (so called global shaking) with 
re-estimation of the neighbor inclusions in the testing window w^ m ^ after each local shaking of the 
inclusion Uj. In so doing, the increasing number of shaking has a twofold effect. On the one hand, the 
system becomes more homogeneous and well-mixed and on the other hand, the stochastic fluctuations of 
statistical parameters (such as e.g. gkmif)) estimated by averaging of these structures tend to diminish. 

The number of points in a sphere of radius r about a randomly chosen point Xj is a random parameter. 
Estimation of its expectation defining K-function leads to loss of statistical information which can be 
presented in more details if one uses a histogram form. Namely, we evaluate the histograms of the 
average number rih of inclusions in the spherical histogram window w£ = {x : |x - x,| < R {k) a^} (i.e. 
the fractions p{nt) of the histogram windows containing the numbers rih of inclusions) as a measure of 
inhomogeneity of an inclusion distribution inside w k . Obviously, the descriptor p(nh) is a more sensitive 
measure of the local statistical homogeneity of the configuration analyzed than the K(r) is. 

Thus, the scheme of the modified CRM (called shaking collective rearrangement model, SCRM) 
accompanied by the shaking procedure is the following. After a few repetitions estimating new velocities 
for the colliding growing inclusions (described in Subsection 3.1), the inclusion shaking (described in 
Subsection 3.2) is carried out. This new algorithm has a few advantages. The SCRM is not as optimal 
as the original CRM for modeling close-packing configurations simply because the added procedure of 
random shaking is just focused on the "destruction" of dense "locked" local configurations in some testing 
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windows that would be manifested from one side by reduction of the first pick of the RDF g(2a + ). From 
other side, shaking leads to the generation of highly homogeneous and mixed structures that can be 
controlled by the histograms of fractions p(n t ) of testing window containing rih particles. Moreover, a 
distinguishing feature of the initial CRM is a large noise of the estimated RDF (at least for the large 
particles, see numerical results in Section 5) because these data should be considered as merely a single 
realization of such a random generation process. In order to provide statistically more reliable results it 
would be necessary either to average several realizations or increase a size of a sample. However, this 
repeating procedure is not necessary in the protocols accompanied by a shaking procedure because a 
configuration generated by a few global shaking can be regarded as a separate realization. 

It should be mentioned that the shaking procedure is a very time consuming one. As it was mentioned, 
the shaking has twofold effect. On one hand side, sequential arrangement of global shaking destructs the 
dense locked local configuration generated by the CRM which is protocol dependent. Optimization of 
this part of algorithm from the point of view of CPU time reduction is hampered. However, because 
of the ideally parallel nature of the problem, the second goal of shaking related with diminishing of 
stochastic fluctuations of statistical parameters estimated can be very effectively reached by using of 
parallel computing at the standard instruction-level parallelism that is performed in the computer codes 
developed. 

5 Numerical results 

5.1 Unimodal distribution 

It should be mentioned that the most experimentally investigated cases of the RDF correspond to 
the close packing (from the loose packs c = 0.59 to dense packs c = 0.64), see van Blaaderen and Wiltzius 
(1995), Gidaspov and Huilin (1998), Aste et al, (2004, 2005), Saadatfar (2009), Kurita and Weeks (2010). 
The volume fractions < c < 0.6 are less well understood, and, moreover, their investigations are usually 
limited by suspensions with the different effects of interparticle interactions (see, e.g., Sierou and Brady, 
2002) that makes difficult their application to the solid composites. Estimation of the RDF from numerical 
simulations of random structures with < c < 0.6 is still less explored (see the use of the hard core model 
for simulation in Segurado and Llorca, 2002). 

The RDF g(r) is well investigated only for identical spherical (3D and 2D cases, see for references, 
e.g., Buryachenko, 2007) inclusions with a radius a. Two alternative analytical RDFs of inclusion will be 
examined for 3D case 

g(r) = H(r-2a) (10) 
g(r) = H(r-2a)[l + gY(c)gW(r)], (11) 

where H denotes the Heaviside step function, r = |x^ — x 9 | is the distance between the nonintersecting 
fixed inclusions Vi and moving one v q , c is the volume fraction of particles of the radius a; g^ (r) = 1 at 
r = 2a + and g^ (c) = g{2a + ) — 1. The formula (10) describes a well-stirred approximation (differing from 
the RDF for a Poisson distribution by the availability of "excluded volume" with the center x^ where 
<?(|xi — x g |) = 0) while Eq. (11), proposed by Willis (1978), takes into account a neighboring order in the 
distribution of the inclusions (see also Stell and Rirvold, 1987): 

gf{c) = g%(c) ee 0.8006(^) , g?{r) = cos (^e 2 ^/-). (12) 

The restriction is that c ^ 1/8, when (c) ^ (12i). For avoiding of this restriction, Kanaun (1990) 
proposed a modification of the formula (12i) when g^ (c) + 1 coincide with the first peak of the popular 
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Percus-Yervik (1958) approximation (see Throop and Bearman, 1965) 

2 + c 



PY i „\ _ 



2(1 -c?) 



-1, 



(13) 



which is known to be most accurate when the density is not too high (c < 0.54, see, e.g., Schaertl and 
Sillescu, 1994). Following Kanaun (1990), we will call Eq. (11) with the representations g^(c) (13) and 
9? ( r ) (1^2) as the Willis approximation (11). 

The analytical RDFs (10), (11) and Percus-Yervik (P-Y, 1958) approximation (see numerical results 
for P-Y approximation in Throop and Bearman, 1965) are compared with the numerical simulations by 
the CRM and SCRM (the number of spheres is n = 2000) for c = 0.575 in Figs. 1 and 2 for the ranges 
of a normalized coordinates 2 < R = r/a < 2.2 and 2.5 < R < 4.5, respectively. As it is well 
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Fig. 1: The RDFs g(r) vs r/a at c = 0.575 estimated Fig. 2. The RDFs g(r) vs r/a at c = 0.575 estimated 
by CRM (1), SCRM (2), Eq. (5.2) (3), P-Y (4) by CRM (1), SCRM (2), Eq. (5.2) (3), P-Y (4) 

known, the values <?,f Y (c) = gY( c ) are underestimations of the corresponding first picks at r = 2a + 
obtained by both the CRM and SCRM. It should be mentioned that the Monte Carlo shaking process is 
both the most structurally influential and, computationally, the most intensive part of the whole packing 
process. So, packing by the CRM for n = 2000 takes a few seconds for a PC with a 2.4 GHz processor 
while one global shaking requires 5 min and stabilization (rather than elimination of a statistical noise of 
descriptors evaluated) of statistical parameter estimations is reached after 10 global shakings. Thus, the 
shaking process is very time-consuming and sensitive to a few competitive variables (such as i?X and 

iV s k ) and the frequency ratio of the alternation of the CRM's iterations and the shaking which should 
be optimized (optimization of the CRM algorithm was considered in Maggi, et ai, 2008). So, increasing 
of the shaking radius leads to both the positive and negative effects. From one side, increasing of 
R^} leads to large amplitude of possible moving of an inclusion i>j but from other side, it brings necessity 
to increase both the near-neighbor list J\f(i) and the limiting shaking number that is very time- 

consuming. The optimal parameters R^ = 0.02 and = 150 were found at the data estimations in 
Fig. 1 and 2. It should be mentioned that derivative in the RDF (4) is estimated by the method of finite 
difference with a chosen Ar. If Ar is too large, g(r) is smoothed and peaks are obscure; if too small 
g(r) has large statistical noise with perhaps spurious peaks because only a small number of particles can 
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be observed in a spherical shell with the thickness Ar and an inner radius r. Some trial runs indicate 
the optimal Ar = 0.01a; decreasing Ar from 0.015a to 0.01a leads to increasing of the first peak at 
r = 2a + on 3% while the different choices of Ar within a broad range of 10 _4 a to 10 _2 a lead to almost 
indistinguishable results. 

It is interesting that the estimations of g(r) (4) can be performed through the finite difference analogs 
of the central (Aste et at, 2004) and right (Maggi, Stafford, and Jackson, 2009) derivatives (Ar > 0) 

g ± (r,Ar) = b[K(r + Ar/2) - K (r - Ar/2)](Ar) -1 , 

5 +(r,Ar) = b[K(r + Ar) - if (r)](Ar)" x , (14) 

respectively, where the constant b = r 1_rf /(dw ( j) is fixed by imposing that asymptotically g(r) -> 1 as 
r — > oo. The formulae (14i) and (142) providing the accuracies 0(Ar 2 ) and O(Ar), respectively, lead to 
close results at A < 10~ 2 a if K (r) is differentiable in r. However, K{r) has discontinuity at r = 2a and 
is not differentiable at r = 2a: K(r — 0) = and K(r + 0) > 0. Therefore, a formal incorrect use of Eq. 
(14i) leads to the dramatic result at r = 2a: 

g ± { r, Ar) = b ^ (? ' + Ar/2) = V(r, Ar/2), (15) 
ZXr 2 

and, therefore at r = 2a 

9 ± (r) = g + (r)/2 (16) 

while <7±(r) = g + { r ) at r > 2a; here g 1 (r) = lim<? (r, Ar) and g + (r) = lim<7 + (r, Ar) for Ar — > 0. 
Apparently the use of Eq. (16) has led Aste et al. (2004) to an unusual power law approximation g(r) oc 
(r — ro)~ Q in the range r < 2.8a with the singularity at ro = 2.06a while a similar behavior, but with 
ro = 2a (that agree with our results), was reported by Silbcrt et al. (2002) for molecular dynamic 
simulations. 

The RDFs estimated from the CRM and SCRM simulations in Figs. 1 and 2 have a well-known form 
(see, e.g., Lochmann, Oger and Stoyan, 2006; Maggi et al., 2008) with the main maxima appearing at 
well-defined positions, namely at r = 2a, 2v / 3a, 4a (the first, second, and third peaks, respectively ) and 
so on (which can be also presented in a nondimcnsional form R — 2, 2\/3, 4) and then g{r) continues to 
fluctuate with decreasing amplitude. The smoothing of the RDFs obtained from the SCRM simulation was 
performed by the averaging over P = 50 realizations of paralleled simulations. These maxima represent 
the distances between the centers of a reference sphere to its first, second, third neighbors etc. The first 
very pronounced peak (discontinuity) at the hard-core distance r = 2a is a result of a direct contact 
of spheres which is a most frequent nearest-neighbor distance. The second peak at r = 2\/3a w 3.464a 
(see Fig. 2, and 3(A)) corresponding to a Face-Centered Cubic (FCC) lattice and appearing at the high 
packing fraction is completely lost at the low volume fraction of particles. The peak at r = 4a (growing 
faster with growing c than the second one) sorts well with the placement of three touching spheres in a 
row (see Fig. 3(B)). The gap distance r = 2.75a = 2.85a corresponding to the empty space between the 
first and second neighbors of typical particle is also clearly observable. The third peak at r = 4a is well 
predicted by the Percus-Yervik (1959) approximation which does not expect the second peak at r = 2\/3a 
for any c. Willis (1978) approximation (11) (curves 3 in Figs. 1 and 2) is in agreement with simulation 
observation (curves 1 and 2) only qualitatively. Arrangement of inclusions presenting in Figs. 3(A) and 
3(B) the fixed left inclusion and defining the second and third, respectively, peaks can be denoted in the 
forms Vi — (pi +i>j) — Vi and Vj — (u,) — Uj, respectively. Here the inclusions placed within the parentheses 
represent the first shell of nearest moving particles. 
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As can be seen in Fig. 2, the SCRM leads to the smoothing of the RDFs g SCRM {r,c). But the 
smoothing of the RDFs (see the curve 1 in Fig. 2) could be reached significantly cheaper in the framework 
of the initial version of the CRM by increasing of either N or the number of realizations X. Nevertheless, 
any smoothed RDF obtained by averaging of RDFs of configurations generated by the CRM don't remove 
the birth-marks of a fundamental shortcoming of the CRM such as protocol dependent and too strong 
local arrangement in the nearest zone (r = 2a + ). 

In doing so, a fundamental advantage of the SCRM is, first of all, a destruction of dense "locked" local 
configurations (which are protocol dependent) leading to the most homogeneous and mixed structures. 
It is especially pronounced in Fig. 1 where SCRM demonstrates a reduction of g(2a + ) estimated by the 
CRM of 43% (destruction of the nearest zone arrangement which is responsible for the local jamming). 
For more detailed analysis of RDF's regularities, one estimates g(r) by both the CRM and SCRM (at the 
same simulation parameters as in Figs. 1 and 2) for more dense packing at c = 0.6 and c = 0.63. As it is 
expected, the CRM's peaks are always greater than the SCRM's ones, significantly so in some cases. For 
example the SCRM reduces the first peaks (r = 2a + ) estimated by the CRM of 43%, 21%, and just 4% 
at c = 0.575, 0.6, and 0.63, respectively. In doing so, the shaking procedure at r > 2.5a leads to nothing- 
more than a smoothing at any volume fractions c of inclusions; so, even for c = 0.6 a systematic reduction 
of the peaks at r = 2\/3a and r = 4a is no better than 10% while there is no systematic difference outside 
of peak's areas at r > 2.5a between RDFs of configurations simulated by both the CRM and SCRM. 
So much significant reduction of g(2a + ) at c = 0.575 and c = 0.60, is explained by the fact that at the 
volume fractions c far from the jamming limit c = c lam = 0.6366 a particle configuration has many local 
areas with relatively small densities where the particles can be removed from the destructed nearest zones 
r = 2a + by the shaking. The last statements is confirmed by Fig. 4 where the histogram of fractions 
p(rih) of histogram window = 2.5a containing inclusions are presented for c = 0.6. As can be seen, 
shaking leads to reduction of the numbers of histogram windows 



r = 2V3a 
< > 




Number of particles 

Fig. 3: Possible sphere configurations and the peak lo- Fig. 4. Histogram of fractions p(n>h) of histogram win- 
cations in the RDF: (A) Four close spheres Vi — (Vi — dow containing nh spheres generated by the CRM (solid 
v i) ~ Vi, (B) Tree spheres in a row Vi — («») — V{ line, 1), SCRM (dotted line, 2) 

with small and largest numbers of particles and increasing of window number with average number of 
spheres. However, at the c > 0.63 the areas with local dense "locked" packing are increasing in the 
volume and occupy in the limiting case c = 0.6366 a full sample. Because of this, shaking has very poor 
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effectiveness at c = 0.63 and reduces the first peak at r = 2a + just by 4%. 

The availability of 3D imaging tools such as computed X-Rays tomography and microtomography, 
confocal optical microscopy or magnetic resonance imaging, combined with the increased power of com- 
puters allows us to process large quantities of data (including the coordinates of particle centers) at 
reasonable costs. It is important to be sure that the packs simulated by proposed virtual methods have 
the same statistics as the real material packs. It is possible, of course, that real and simulated packs will 
have different statistics according to the manner in which they are manufactured and generated, a matter 
that we are not in a position to address at the present time. But Aste and his colleagues at the Australian 
National University in Canberra have measured the RDFs for a number of very large monodisperse packs 
(~ 150,000 spheres) using X-ray tomography (Aste et ai, 2004, 2005), and have provided us with their 
raw data so that we might make comparisons, and we do that here. Our packs use 35,000 spheres in a 
periodic cube and are generated using the T algorithm (see for details Maggi et ai, 2008) with packing 
fractions of 59.25%; the experimental packing fractions are 59.3%. The comparisons are shown in Figs. 5 
and 6. As can be seen, the curves simulated by the SCRM (n = 2000; P = 80) is better 
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Fig. 5: The RDFs g(r) vs r/a at c = 0.593 estimated Fig. 6. The RDFs g(r) vs r/a at c = 0.593 estimated 
by CRM (1), SCRM (2), Aste's data (3) by CRM (1), SCRM (2), Aste's data (3) 

correlated with experimental data by Aste et al. (2004, 2005) (see also Silbert et ai, 2002) than CRM's 
simulation. In doing so, there is no reason for expectation that RDFs estimated from experimental data 
are the unbiased ones which can be correctly compared with the RDFs evaluated from the simulated data. 
Indeed, both the CRM's and SCRM's data were simulated for unimodal distribution with identical sphere 
size a while for real polydispcrse, non-perfect, spheroidal particles, the distance between the centers of 
particles in contact is not a fixed value 2a but instead it is distributed around an average value. However, 
exploiting of Eqs. (2) and (4) implies the use of some fixed radius <2f of particles which was estimated from 
experimentally found volume fraction c and the intensity n according to the formula 47ra 3 n/3 = c, or, by 
other words, af = (a 3 (C)) 1 ^ 3 , where a(£) is a random radius of particles. This defines a bias of the RDF 
evaluated from experimental data. So, the use of a f = 1.05(a 3 (C)) 1/3 and a f = 0.95(a 3 (C)) 1/3 instead of 
af = (a 3 (C)) 1 ^ 3 leads to the variation of the first peak g(2a + ) on 40% (reduction) and 7% (increasing), 
respectively. 

It should be mentioned that for the SCRM's configurations, the particle reorganization induced 
by shaking is subjected only to geometrical constraints, whereas for real structures the packing is far 
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more complicated and controlled by elastic, hydrodynamic, cohesive forces, and chemical specificity etc. 
(see, e.g., Alcaraz et al., 2008). Our simulation technique is able to isolate the fundamental geometrical 
constraints from other physio-mechanical and chemical effects and therefore the results provide a valuable 
benchmark for evaluating sophisticated packing schemes used in the modeling of real composite materials. 
It is shown that a sufficiently intensive shaking process leads to the stabilization of a statistical distribution 
of the simulated structure that is most homogeneous, highly mixed and protocol independent (in sense 
that the statistical parameters estimated do not depend on either the basic simulated algorithm such as, 
e.g., CRM or manufacturing parameters of the real material production). It is expected (and proved for 
2D case by Buryachenko et at, 2003) that shaking procedure will reduce the peaks (having a clustered 
nature) of experimentally estimated RDF's (by the use of, e.g., X-ray tomography, see Aste et ai, 2004, 
2005) rather than only the peaks of RDFs generated by the CRM (the arguments justified the last 
statement are plausible rather than rigorous and require additional investigations). This inconsistency 
between the real and SCRM's structures is not a deficiency of the SCRM because a subsequent use of 
any advanced micromechanical model (see for references Buryachenko, 2007) assumes a stationary of a 
process X rather than its slightly clustered configuration as for real materials. 

Due to a costliness of the SCRM, it can be useful to extrapolate a value g SCRM ( r ) estimated 
at the volume fraction c\ to other particle concentration c 2 . It can be easy done if one assumes that 
<7 M (r) has in vicinity c = d the same functional dependence on c and r as Willis's approximation 
'(12): g SCRM (r,c) = H(r-2a)[l+g s c CRM (c)-gf CRM ' (r)\ but with the different g s c CRM (c) and g^ CRM (r). 
Assuming also that g s c c ' RM '(d) / ' gf c mM '(c 2 ) = (ci) / g™ (c 2 ) , we can obtain 

9 SCRM (r,c 2 ) - 1 = [g SCRM (r,d) 1] ■ [g?(c 2 ) l]/b c w (d) - 1]. (17) 



It should be mentioned that the RDF was only defined by Eq. (4) for statistically homogeneous 
and isotropic fields of inclusion centers Xi,X2, . . .. Then / 2 (x — y) and g(r) are the continuous isotropic 
functions. However, for the lattice structures, / 2 (x — y) is a sum of a delta functions (with the support at 
the particle centers) and correctness of using of Eq. (4) is not obvious because in such a case estimators 
of g(r) for a stationary ergodic process are difficult to interpret (see Picka, 2007). However, due to the 
practical importance of lattice structure investigations, we present some correct interpretations of Eq. 
(4) for these cases. Namely, we consider as an example a simple cubic lattice (t — x,y,z): t% = i X 2a 
(i = 0,±1,±2, ...,±ni, 2a = l/(m + 0.5). At first we estimated g + (r,Ar) (14 2 ) for ni = 1,2,5 and 
the normalized radius R = r/a and the radius steps Ai? = Ar/a = 0.1,0.01,0.005,0.003: we get a 
picture similar to Nolan and Kavanagh (1995) with g + (2,AR) — 9.55,95.5,190.9, and 318, respectively, 
i.e. g + (2, AR) — > oo for Ai? — > that is intrinsic for the delta function, e.g. S(r — 2a) = dH(r — 2a)/dr. 
However, any integral over the domain containing R — 2 eliminates the delta function that gives different 
opportunities for interpretation of g + (r). In particular, the integral over the shell with the internal radius 
r and the thickness Ar 

yr+Ar 

I(r,Ar)=4irn r 2 g + (r, Ar)dr = [Annr 2 Ar] [g + (r, Ar)}, (18) 

J r 

(where n — N/w) is a number of inclusion centers in the shell. It was detected that for any m (at least 
for m = 1, 2, 5) and any Ai? (at least for Ai? = 0.1, 0.01, 0.005, 0.003), the integers I(r, Ar) = I(r) does 
not depend on n; and Ar. Thus, both items in the square brackets in Eq. (18) depend on Ar but its 
product i(r, Ar) does not depend on Ar (at least at Ar < 0.1). A plot I(r) vs r/a is presented in Fig. 7. 
A visually continuous curve I(r) vs r/a (instead of discontinuous integer function in analytical solution) 
is obtained due to numerical discretization of the delta-function. As we can see, the first peak coincides 
with the number of the nearest neighbors in the simple cubic packing being considered. Thus, the integral 
i(r) uniquely defined has an obvious physical meaning (in contrast to g(r)). 
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Fig. 7. I(r) vs r/a for a simple cubic lattice 



5.3 Bimodal distribution 



For polydisperse structures the counterpart to the RDF is less investigated than for monodisperse 
case and described by the partial RDF gij(r) which is the normalized (by the condition 5^(00) = 1) 
number density of inclusions v^) of type j that are within shells r to r + Ar centered on the inclusions 
uW of type i. gij(r) were estimated by the modified CRM by Maggi et al. (2008) for the size ratio A of 
spherical inclusions more than 6.5. As it was expected, for the large inclusions, 511(7") is stochastically 
very noisy, although trends are apparent. Because of this, as we have already noted, incorporation of the 
shaking procedure in the proposed algorithm is very effective for smoothing of 5^ (r) . 

It is clear that we have to measure the distances between the particles Vi C v^ k ' (ctj = offi) and 
Vj C (a,j — a*-™- 1 ) as the excess distances with respect to contact, that is Ar con = r — {a% + %■); 

rather than the distances r themselves. Then the term H(r — 2a) in Eqs. (10) and (11) should be 
replaced by H(r — a, — aj) Then, we have essentially two ways of scaling these excess distances. One 
is to take the same scaling parameter for the distances in all partial RDFs g^ m (k,m = 1,2) of the 
mixture. This is the situation in the decoupling approximation (see Kotlarchyk and Chen, 1983; Barrio 
and Solana, 2005; see also Lochmann, Oger and Stoyan, 2006) in which the scaling parameter is unity, 
but we can take equally use other suitable scaling parameter such as a^ m K In such a case a different 
scaling distance of each partial RDF for the second peak generated by the second shell of moving particles 
Vi C i>( m ) are placed approximately at the same positions, although the height of the peaks corresponding 
to different RDFs are different. This is particularly remarkable for the first peak at r = af + at. Thus, 
we will present the partial RDFs g^ m as the functions of a scaling parameter Rkm = (r — a 
In a general case, the peaks are a consequence of the bidisperse morphology, but some of them can 
be explained within the monodisperse framework. For example, the coordinates of the first peak of 
gkm and g m k for k 7^ ra coincide in both arguments Rkm = 2 and r while the second 

and third peaks depends on additional information. The nature of the third peak at Rkm = 4 (or 
) described within the monodisperse framework reflects the placement of the touching 
inclusions Vi — (vj) — Vj («, C v^ k \ Vj C t/™)) in a raw (see Fig. 3B). However, a specific bidisperse 
morphology can be manifested by a lot of additional specific configurations. So, it can be generated by the 
configuration m — (wj) — Vj of touching inclusions in a raw that produces a new peak at r = 3a^ + a^" 1 ' 
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or, in normalized coordinates, at R km = 2(1 + Afc m ) (Afc m = /a^). Analogously, the related peaks of 
g m k(r) corresponding to the configurations Vj — (vt) — and Vj — (vj) — vi are placed at r = + 3a^ m ^ 
and r = 3a^ m ^ + a^ k \ respectively, that matches with the peak coordinates of gkm(r)- Thus, the third 
peak (R = 4) of unimodal pack has two possible counterparts of their bidisperse pack with Rkm = 4 and 
Rkm = 2(1 + Xkm) coinciding only at Xkm = 1. In a similar manner, the second peak of unimodal pack 
with a symbolic configuration Vi — (vi+Vi) — Vi (see Fig. 3A) have three counterparts for the fixed inclusion 
V{ C v^: Vi — (vi + Vi) — Vi, Vi — (vj + Vj) — Vj, and Vi — (vi + vj) — Vj. For example, the configurations 
Vi— (vi+Vi)—Vj &n.dvi — (vj+Vj) — Vj corresponding to the peaks Rkm = (v^— l)Afc m + \A + 2Afc m + l and 
Rkm = v3 — Xkm + \J X\ rn + 2Xkm + 1 , respectively. Of course for other the volume fractions c, size ratio 
Xk m , and volume fraction ratio £km = c(*'/c( m ' (k ^ to) being considered, one obtains different results 
for the peak locations (see, e.g., Lochmann, Oger and Stoyan, 2006) when the different "monodisperse" 
and "bidisperse" peaks are pronounced in varying degrees. So, visualization of peaks in unimodal packs is 
only managed by a single parameter c while the peaks of bimodal packs are managed by three parameters 
c, Xkm, and £ fem (fc ^ to). 

Due to the absence of the analytical representations of gij(r) = 1,2) for polydispcrse structures, 
it is possible to construct <?y(r) within monodisperse considered framework (10) and (11) in the following 
manner {Rij = (r — ai)/aj + 1): 

g ij (r)=H(R ij -2), (19) 
9ij (r) = H{R. l0 - 2) [l + c) - 1] cos( 7 rii ii )e 2 ( 2 - fl «)] , (20) 

where c = c^ 1 ' + c^ 2 '. The RDFs contact values fl£j(£, c) were obtained by Attard (1989) from Ornstein- 
Zernike equation using Pcrcus-Yevick closure of the pair and triplet levels. In particular, the P-Y result 
obtained for the bulk mixture at the pair level is 

ffijKj) = [ajduiaid + CLigjjiaj^/aij, (21) 
Qiiiau) = [l-C 3 + 3a 4 C 2 ](l-C 3 ) 2 , (22) 

where = a; + aj is a mutual diameter, and £ a = (7r2 Q_1 /3) J2i a ? ■ Taking the limit as c' 2 ' —> 0, so 
that £3 = one has 

= i 1 + c( ) 7A^TTyJ (l-c^))" (23) 

which reproduces the respective limit 512(012) = (1 + c /2)/(l — c^ 1 ') -2 corresponding to the usual P-Y 
value for unimodal pack g^ Y (c) (20). 

The paper by Chong, Christiansen and Baer (1971) is dedicated to the experimental measurement 
of effective Newtonian suspensions rj*/r]^ ' of spherical bimodal particles with the size ratio A = A21 = 
a (2)/ a (i) = 0.477,0.313, 0.138, and the volume fraction ratio £ = ^ 2 i = c( 2 '/c (1) = 0.333. Because of 
this one considers the partial RDF's gij{Rij) (no summation over i,j — 1,2) simulation for bimodal 
packs with A = 0.313 and c = 0.6. The typical first peak at Rijh — 2 of the gij(Rij) presented for 
the RDF gii(i?n) in Fig. 8 is analogous to the first peak for the unimodal pack (see Fig. 1); shaking 
parallel procedure SCRM yields a reduction of the first CRM's peak on 25%. The curves gij(Rij) (see 
the curves 2 in Figs. 9-11) behave at Rij > 2.2 less trivial than g(R) (see Fig. 2). So the second peak 
-R111 2 = 2-64 holds a clearly defined bidisperse nature and corresponds to the configuration V\ — (1*2) — v i- 
R\i\2 = 2.64 ~ 1 + 2A + 1. The "monodisperse" peaks Rn\z = 3.464 and Rnu = 4.0 correspond to the 
peaks R = 2\/3 ~ 3.462 and R = 4 of the unimodal pack in Fig 2. The pronounced peak in Fig. 9 at 
i?ii|s = 4.54 corresponds to the case when the shell either — (u 2 + V2)— or — (ua)— is wedged in between 
the three spheres in a row v% — («i) — V\ that moves the peak i?iiu = 4.0 on either 2(^1 + 2A — 1) ~ 0.55 
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Fig. 8: E*(c)/E lS)) vs c estimated by the MEFM with 
RDF simulated by CRM (1), SCRM (2), Willis approx- 
imation (3), well-stirred approximation (4), and by the 
MTM (5); o - experimental data by Smith (1976) 



Fig. 9. /i*(c)///°) vs c estimated by the MEFM with 
RDF simulated by Willis approximation (1), SCRM's 
approximation (5.5) (2), CRM (3), SCRM (4), PY (5), 
well-stirred approximation (6), and by the MTM (7); o 
- experimental data by Kreger (1972) 



first and the second neighbors of typical particle v\ and symbolically denoted by the form v\ — (vi +v\) ^ 
at R — 2 + y3/2 — 2.87. This gap is manifested in g(R) but not in 912(^12) an d 322(^22) where the first 
gap distance is formed by a bidisperse configuration v\ ~ (vi) and V2 — (t>i) ~— ' at R12 = 1/A ~ 3.2 
and R22 =1/A^ 3.2 (observed as 3.3 in Figs. 10 and 11), respectively, rather than by a configuration 
Vi — (V2 + V2) and v 2 — (V2 + V2) at R12, R22 = 2.87. The distance i?i2, -R22 = 1/A ~ 3.2 corresponds 
to the middle surface of the shell — (i>i)— surrounding a typical fixed point v± (or V2), occupied by the 
centers of moving particles V\ and inaccessible to the centers i>2- Moreover, a similar gap at the distance 
Rn = 2 + A ~ 2.3 (see Fig. 9) of gu(Rn) holds the same bidisperse nature and is produced by the 
shell —(V2)— surrounding a typical point V\ and unattainable to the centers of moving particles V\. The 
next (after i? 12 = 4) most pronounced peak of 312(^12) in Fig. 10 is analogous to the peak i? 12 = 4 
with the configuration v\ — (1)2) — «2 where a middle sphere — (1*2) — is replaced by a sphere — (i>i)— in 
the configuration V\ — (fi) — V2 that leads to displacement of the peak R\2 = 4 on 2/ A — 2 to the right 
to R12 = 2 + 2/A ~ 8.3. In a similar manner, the peak of 322(^22) in Fig. 11 at i?22 = 2 + 2/A ~ 8.3 
produced by the configuration V2 — (vi) — V2 is clearly defined. 

The RDFs curves gu(Rn), (712(^12) and 322(^22) in Figs. 8-11 obtained for N = 2000 were smoothed 
by averaging over P = 80 realizations of the paralleled SCRM's generations. For analysis of the impact of 
a sample size N, one performs a similar modeling for ./V = 7000 for M = 3 global shaking and compares 
it with analogous simulations for N — 2000 with M = 3 and M — 60 global shaking (see Figs. 12 and 13). 
As can be seen the RDFs estimated for the samples with N = 2000, M = 60 and N = 7000, M = 3 are 
close while the difference of the RDFs estimated for TV = 2000 with M = 3 and M = 60 global shaking is 
significant. The results of comparison indicate on a dual nature of smoothing averaging which accuracy 
can be increased by two ways either by increasing of the number of realizations (the number of paralleled 
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estimations P) or by a growing of a sample size (the number TV). It means that a sample with N = 2000 
can be considered as a statistically representative of an infinite sample. A remarkable difference between 
the RDFs estimated for N = 2000 with M = 3 and M = 60 global shaking points to the fact that M = 3 
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global shaking is not enough for destruction of locked local arrangements generated by the CRM. So, the 
difference of gn(2 + ) for M = 3 and M = 60 equals 12% while a similar dissimilarity for M = 30 and 
M = 60 does not exceed 0.2%. Thus, the results of RDFs estimations are statistically representative in 
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both sense the sample size {N = 2000), statistical stabilization of modeling results (M = 60), and its 
smoothness (P = 80) although the number of big particles is moderately small (N^ = 168 and 
ATM = 587 for N = 2000 and N = 7000, respectively). 

6 Conclusion 

Thus the SCRM proposed makes it possible to generate bimodal random structures of spheres for 
a whole range of particle volume fractions which are primarily homogeneous, highly mixed and protocol 
independent (in the sense that statistical parameters estimated do not depend on the basic simulated 
algorithm such as, CRM). Generalization of the proposed algorithm to any number of fractions is straight- 
forward. 

It should be mentioned that one of the basic challenges of the rocket industry is the optimization of 
size ratios of packs which are usually performed experimentally. Only at the present time this problem can 
be solved numerically. The problem is reduced to a minimization problem f{c^\ . . . ,c^ N ') — > min(max) 
at the constraint cW + . . . + & N ' = c = c max . The objective function f(c^\ . . . , c^ N ') can be described 
as, e.g., either a jamming limit, the values of the RDFs gij(a^ +a^') — 1, . . . ,N) estimated by the 
parallelized SCRM or the effective properties and stress concentrator factors of local stresses estimated in 
an accompanying paper by Buryachenko (2012). However, more detailed consideration of this optimization 
problem is beyond the scope of the current paper and will be analyzed in forthcoming publications by 
the authors. 
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